Exchange-correlation kernel for perturbation dependent auxiliary functions in auxiliary density perturbation theory

Context Analytic exchange-correlation kernel formulations are of the outermost importance for density functional theory (DFT) perturbation calculations. In this paper, the working equation for the exchange-correlation kernel of the generalized gradient approximation (GGA) for perturbation dependent auxiliary functions is derived and discussed in the framework of auxiliary density functional theory (ADFT). The presented new formulation is extended to the unrestricted approach, too. A comprehensive discussion of the implementation of the GGA ADFT kernel, using either the native exchange-correlation functional implementations in deMon2k or the ones from the LibXC library, is given. Calculations with analytic exchange-correlation kernels are compared to their finite difference counterparts. The obtained results are in quantitative agreement. Nevertheless, analytic GGA ADFT kernel implementations show substantial improvement in the computational performance. Similar results are reported for analytic second derivatives of effective core potential (ECP) and model core potential (MCP) matrix elements when compared to their finite difference counterparts in molecular frequency analyses. Method All calculations are performed in the framework of ADFT as implemented in deMon2k. In the ADFT analytic frequency calculations, auxiliary density perturbation theory was used. The underlying two-center exchange-correlation kernel matrix elements are calculated by numerical integration either with analytic or finite difference kernel expressions. Validation calculations are performed with the VWN and PBE functionals employing DFT-optimized DZVP basis sets in conjunction with automatically generated GEN-A2 auxiliary density function sets. In the (Pt3Cu)n cluster benchmark calculations, the RPBE functional was used. For Pt atoms, the quasi-relativistic LANL2DZ effective core potential with the corresponding valence basis set was employed, whereas for Cu atoms, the all-electron DFT-optimized TZVP basis was applied. The auxiliary density was expanded by the automatically generated GEN-A2* auxiliary function set. We run all benchmark calculations in parallel on 24 cores.


Introduction
Over the last three decades, Kohn-Sham density functional theory (DFT) [1,2] methods have become the workhorse in computational chemistry and materials science [3].Their excellent performance-to-accuracy ratio permits reliable energy and structure predictions for systems with several hundreds of atoms.Although there exists still an accuracy Dedicated to Professor Toro-Labbé on the occasion of his 70th birthday.
Extended author information available on the last page of the article gap to the targeted chemical accuracy of 1 kcal/mol, it is foreseeable that this gap will close with the further development of Kohn-Sham DFT.To this end, it is important to improve density functional approximations (DFAs) such that the overall computational performance is not jeopardized.A possible framework for such developments is auxiliary density functional theory (ADFT) [4] which is based on the variational density fitting for Coulomb [5,6] and Fock [7,8] energies and uses the resulting auxiliary density for the evaluation of the exchange-correlation contributions [9,10].As a result, ADFT local density approximation (LDA) and generalized gradient approximation (GGA) as well as hybrid calculations have a formal cubic scaling with respect to the number of basis functions.No four-center electron repulsion integral (ERI) evaluations nor basis function product evaluations at grid points are needed in ADFT.Despite this enormous computational simplification, ADFT provides for a given DFA the same accuracy in terms of structure parameters, relative energies, and harmonic frequencies, to name a few, as corresponding Kohn-Sham calculations.Therefore, ADFT is particularly well suited for the reliable optimization of complex nanometric systems like large transition metal clusters [11,12] or Born-Oppenheimer molecular dynamics simulations [13,14].
Because the ADFT energy expression is variational, analytic energy derivatives can be straightforwardly calculated, e.g., energy gradients by the 2n+1 theorem of perturbation theory [15].This permits an alternative perturbation theory formulation to the commonly employed coupled-perturbed Kohn-Sham (CPKS) approach [16].We have named this alternative approach auxiliary density perturbation theory (ADPT) [17,18] to emphasize its close connection to ADFT.Since its initial formulation, ADPT has been used to calculate polarizabilities [19][20][21][22][23], Fukui functions [24,25], electron binding energies [26], alchemical derivatives [27], and various magnetic properties [28][29][30][31][32][33][34].In all cases, CPKS equivalent results were obtained with a much reduced computational demand.In fact, for LDA and GGA perturbation calculations, timings comparable to the ADFT selfconsistent field (SCF) approach are observed.More recently, this has been extended to time-dependent ADFT excited state calculations obtaining similar computational performance [35][36][37].One reason for the improved computational performance of ADPT arises from the numerical exchangecorrelation kernel calculations, i.e., the second functional derivatives of the exchange-correlation energy.In ADPT, these matrix elements contain only auxiliary functions from two centers.As a result, the corresponding kernel matrix has only the dimension of the number of auxiliary functions.Key ingredients for the efficient exchange-correlation kernel calculations are analytic formulas for the various LDA and GGA functionals.Whereas the LDA ADFT kernel formulas are rather straightforward to derive, the corresponding GGA formulas are more involved and deserve some attention.For perturbation-independent auxiliary functions, the generic formulas for closed-shell ADFT kernel calculations have been published a few years ago [38].However, the recent development of ADFT second analytic energy derivatives [39] employing ADPT requires corresponding kernel calculations for perturbation dependent auxiliary functions.Again, computationally efficient formulations are needed because we aim for frequency analyses of systems with many hundreds to several thousands of atoms.Such calculations are not only mandatory for the characterization of the optimized structures but also for initial Hessian matrix calculations, either for structure minimizations or transition state optimizations.Again, the focus is on GGA kernel formulas which we will derive in this work.To keep the discussion most general, we present the GGA kernel working equations in the framework of the unrestricted formulation.Furthermore, we also validate the recently implemented LibXC interface [40] of deMon2k in the context of these kernel calculations.Because we found in some transition metal cluster frequency analyses a computational bottleneck arising from the finite difference calculations of integrals for effective core potential (ECP) and model core potential (MCP) second derivatives, we also report the corresponding analytic second derivative formulas that were newly implemented into deMon2k [41] in the framework of this study.
The paper is organized as follows.The subsequent theory section contains four subsections.After a brief introduction of second analytic ADFT energy derivatives, the generic working formulas for the GGA potential and kernel calculations for perturbation dependent auxiliary functions are presented.The two following subsections discuss the finite difference GGA ADFT kernel formulations that can be employed if explicit kernel expressions for a GGA functional are missing, and the analytical second derivative formulas for ECP and MCP matrix elements.Section 3 provides the computational details for the validation and benchmark calculations that we present and discuss in the following section.All calculations are performed with the deMon2k code using either the native or LibXC exchange-correlation functional implementations.Finally, in Section 5, the conclusions are summarized.

Theory
In the linear combination of Gaussian-type orbitals (LCGTO) approximation, neglecting for the sake of simplicity explicit spin dependency, the second analytic ADFT energy derivatives with respect to atomic coordinates λ and η are given by the following [39]: In Eq. ( 1), superscripts in parentheses denote derivatives with respect to the corresponding atomic coordinates λ and η.Thus, E (λη) corresponds to a Hessian matrix element.The Greek letters μ and ν denote (contracted) atomic GTOs, whereas the Latin letters with a bar, k and l, represent primitive atom-centered Hermite Gaussian auxiliary functions.Furthermore, P μν is an element of the (closedshell) density matrix, and H μν is an element of the core matrix, incorporating kinetic, nuclear attraction, and any other external potential contributions, e.g., external electric or magnetic fields.The symbol || in the three-center ERI shorthand notation denotes the Coulomb operator 1/|r 1 −r 2 |.It also separates functions from electron 1, in the bra, from those of electron 2, in the ket.The Coulomb and exchangecorrelation fitting coefficients are denoted by x k and z k , respectively.The matrix elements S μν , W μν and G k l belong to the overlap, energy-weighted density, and Coulomb matrices.Superscripts on matrix elements or fitting coefficients indicate derivatives with respect to corresponding atomic coordinates.

Analytic GGA exchange-correlation potential
Of particular interest to our further discussion are the matrix elements in Eq. (1) that contain v xc and v (η) xc , which are the auxiliary density exchange-correlation potential and its derivatives that result in the corresponding kernel expressions.Note that in ADFT, all these expressions are evaluated with the auxiliary density obtained from the variational fitting of Coulomb and Fock energies.In the case of a GGA functional, the exchange-correlation energy takes the following form: In Eq. ( 2), e xc ( ρ, γ ) denotes the density-weighted exchange-correlation energy density depending on the auxiliary density, ρ(r), and its gradient square, γ (r).For simplicity of notation, we omit the explicit position dependencies from the exchange-correlation energy density, potential, and kernel throughout the discussion.In the unrestricted approach, the auxiliary density is expanded as follows: The x α k and x β k are the spin-polarized Coulomb fitting coefficients, obtained from separate spin-dependent fitting equations.The density gradient corrections are included by the scalar γ σ τ (r) field, with σ and τ being labels for the spin, either α or β: Therefore, the explicit form of the unrestricted exchangecorrelation energy in Eq. ( 2) is given by the following: The first derivative of the unrestricted exchange-correlation energy, Eq. ( 5), with respect to the atomic coordinate λ yields the following: For the sake of clarity of presentation, we use explicit summation in Eq. (6).The shorthand notation ρσ u (r) stands for ∂ ρσ /∂u with u being x, y, and z.Thus, ρσ u (r) denotes the components of the auxiliary density gradient, ∇ ρσ (r), with respect to the electronic coordinates.Expanding the auxiliary density and density gradient derivatives as, and yields the following: Defining the unrestricted GGA ADFT exchange-correlation potential operator [38], allows the following shorthand notation for the exchangecorrelation energy derivative of Eq. ( 6): Note that the first term of Eq. ( 11) will be absorbed in the Pulay term of the ADFT energy gradients [42] and, therefore, will not be explicitly calculated.Terms analogous to the second sum of Eq. ( 11) appear in the Hessian matrix elements given in Eq. ( 1) in the form of the k(λ) |v xc [ ρ] and k(λη) |v xc [ ρ] matrix elements.Whereas the first terms are explicitly evaluated in the response ADPT calculations, the second term's contributions to the skeleton Hessian matrix are indirectly evaluated by translational invariance.This permits the use of moderate grids for the numerical integration in the Hessian matrix calculation even without the explicit inclusion of atomic weight function derivatives [43].
For the implementation in deMon2k, it is more convenient to formulate the exchange-correlation potential of Eq. (10) in terms of ρσ (r) and γ σ τ (r).To this end, we apply the chain rule to Eq. ( 10) obtaining the following: For the derivative of γ σ τ (r) with respect to ρσ u (r) holds: Inserting Eq. ( 13) into Eq.( 12) and the result into Eq.( 11) yields the explicit form of the unrestricted exchangecorrelation energy derivative as implemented in deMon2k: The corresponding exchange-correlation potential, v σ xc [ ρ, γ ], is given by the following:

Analytic GGA exchange-correlation kernel
We now turn to the analytic GGA ADFT exchange-correlation kernel calculation in ADPT.To this end, we derive Eq. ( 11) with respect to a second atomic coordinate keeping in mind that the first term of Eq. ( 11) is absorbed by the Pulay relation.Thus, only the second sum of Eq. ( 11) needs to be explicitly considered: The first two terms of Eq. ( 16) are already discussed in Section 2.1 and can be calculated by the exchange-correlation potential expressions given above.Thus, we focus in the following only on the third term of Eq. ( 16).
To proceed, we derive the exchange-correlation potential with respect to the atomic coordinate η: Expanding the auxiliary density and density gradient derivatives according to Eqs. ( 7) and ( 8) results into the following: Inserting the explicit expression of the GGA exchangecorrelation potential from Eq. ( 10) into Eq.( 18) yields the following: At this point, it is convenient to define the GGA ADFT exchange-correlation kernel according to Eq. ( 19).Therefore, we obtain as GGA kernel expression in ADFT: As Eq. ( 20) shows, the exchange-correlation kernel has in general four spin contributions, namely αα, αβ, βα, and ββ.As developed in Eq. ( 20), the kernel is also non-local, i.e., depends from r and r .However, for LDA and GGA, e xc ( ρ, γ ) is an ordinary function.As a result, the variables r and r collapse in the integration over the exchangecorrelation kernel as outlined in [38] and [44].
Back-substitution of the kernel definition, Eq. ( 20), into the exchange-correlation potential derivative, Eq. ( 19), yields for the matrix element in the last term of Eq. ( 16): Therefore, we find explicit contributions to the ADFT Hessian matrix elements from the second derivatives of the exchange-correlation energy functional with respect to the atomic coordinates η and λ: For the kernel implementation in deMon2k, it is more convenient to formulate the exchange-correlation kernel of Eq. ( 20) in terms of ρσ (r) and γ σ τ (r).Applying the chain rule to the terms involving the σ and τ density gradients in Eq. ( 19) and using Eq. ( 13), we obtain the following: Hence, the implemented GGA kernel formulation in ADFT is given by the following:

Finite difference GGA exchange-correlation kernel
In the last section, the working formula for the unrestricted GGA ADFT kernel, Eq. ( 24), in deMon2k was presented.This implementation requires second derivatives of the density-weighted exchange-correlation energy density.
In cases where these derivatives are not available, the matrix elements of the GGA ADFT kernel can be calculated by a finite difference formulation.This approach requires only the implementation of the corresponding GGA exchangecorrelation potential formulas that are mandatory for SCF calculations.To this end, the exchange-correlation potential is evaluated at symmetric arbitrarily small changes in the auxiliary density at a given grid point.Following the definition of derivatives, the unrestricted finite difference formulation of the GGA ADFT kernel can be expressed as follows: Equation ( 25) represents the arbitrary changes in the auxiliary density by ± k(λ) (r), where is the step size, and k(λ) (r) is the derivative of the auxiliary function k with respect to the atomic coordinate λ.This choice is made with the intention of constructing the kernel integrals appearing in Eq. ( 22).
Take as an example the kernel integral in the third term of Eq. ( 22).It can be expressed by the following finite difference formulation: Equation ( 26) can be used to calculate the exchangecorrelation kernel integrals appearing in Eq. ( 22) for any GGA functional for which exchange-correlation potential formulas are implemented.The integral in Eq. ( 26) is evaluated numerically with the same grid used in the SCF.In deMon2k, the step size is set to 10 -9 a.u., which yields consistent and robust results [20].

Analytic second derivatives of ECPs and MCPs
In the original frequency analysis module of deMon2k, the second derivatives of ECPs and MCPs are calculated by finite differences from the corresponding analytic gradients.For systems with many ECP or MCP centers, these finite difference calculations can become a computational bottleneck.Therefore, we implemented within this work analytic second derivatives for ECPs and MCPs.The relevant contributions to the Hessian matrix from second derivatives of ECPs or MCPs, U , centered on atom C have the general form: Taking r C = r − C for ECPs, the pseudopotential operator is given by the following: The first term on the right-hand side is a purely local central force potential located on atom C. The second term is a double sum of semi-local operators with a fully local radial part and non-local angular part, here expressed by normalized real spherical harmonics, S lm (r C ).In this notation, rC is the unit vector in the direction of r C .For MCPs, a similar operator form is used: In Eq. ( 29), B l is a constant and φ lm (r C ) is a core orbital.Since core orbitals are normalized, the terms inside the sum are true projectors, fully non-local.Analytical calculation of second derivatives was implemented taking advantage of translational invariance.Differentiation of the operators, especially the semi-local ECP operator, is not advisable.Therefore, explicit differentiation of the pseudopotential operators was avoided, only basis functions are differentiated.
The Kronecker delta in Eq. ( 30) tests only if λ or η correspond to the same atom where the differentiated atomic basis function is located.The derivative of a Cartesian Gaussian function is a combination of Gaussians with shifted angularity: ± 1 for first derivatives and ± 2 and 0 for second derivatives.Avoiding to differentiate the pseudopotential operator allows us to apply directly the semi-numerical integration algorithm developed for ECPs [45].Evaluation of semi-local ECPs and local ECPs and MCPs' second derivatives has a computational cost slightly larger than the evaluation of the corresponding energy integrals.The extra cost arises from the increased angularity of the required angular integrals which are, however, independent of the number of atoms.Non-local projectors in MCPs are easily decomposed into core-valence overlap integrals, and their derivatives can be evaluated very efficiently since they are completely analog to the overlap matrix derivatives.

Computational details
The open-shell analytical GGA exchange-correlation kernel formula, as detailed in the preceding section, Eq. ( 24), has been incorporated into the LCGTO-ADFT code deMon2k [46].For the kernel calculations, either the native or LibXC partial derivatives of e xc ( ρ, γ ) are used.We denote them by native or LibXC in the following.The test calculations are performed with the Vosko, Wilk, and Nusair (VWN) LDA functional [47] in conjunction with Dirac exchange [48] or the GGA Perdew-Burke-Ernzerhof (PBE) functional [49] as implemented in deMon2k or LibXC.The Kohn-Sham orbitals are expanded with the double-zeta valence polarized (DZVP) basis set [50], while the auxiliary density is computed using the automatically generated GEN-A2 auxiliary function set [51].Thus, the variational fitting of the Coulomb potential [5] is used, and the exchange-correlation energies, potentials, and kernels are calculated with the resulting auxiliary density.This ADFT approach is also used for the analytic frequency calculations.In all calculations, the self-consistent field (SCF) energy and auxiliary density convergence criteria are set to 10 -5 a.u. and 5×10 -4 a.u., respectively.For the numerical integration of the exchange-correlation energy and its derivatives, an adaptive grid with a grid tolerance of 10 -5 a.u. was used.
As benchmark calculations, we performed harmonic frequency analyses of (Pt 3 Cu) n clusters with n being 1 to 11.The calculations are executed with the RPBE [52] GGA functional, again using either the native or LibXC implementation of it.In these calculations, the extended GEN-A2* auxiliary function sets [53] are used.For the platinum atoms, a doublezeta basis set from the Los Alamos National Laboratory was employed in conjunction with an 18-electron quasirelativistic effective core potential (QECP| LANL2DZ) [54].For the copper atoms, an all-electron triple-zeta valence polarized (TZVP) basis set optimized for GGA functionals [51] was used.Initial structures of the (Pt 3 Cu) n clusters were taken from [11,12] and tightly optimized using an SCF energy convergence criteria and an auxiliary density convergence criteria of 5×10 -8 a.u. and 10 -5 a.u., respectively.The structure optimization tolerance was set to 3×10 -4 a.u.for the root-mean-square forces.The frequencies are calculated using the optimized cluster structures by analytic second ADFT energy derivatives employing ADPT.The required exchange-correlation kernels were calculated with the analytical and finite difference methods.In these calculations, a fixed fine grid was used for the numerical integration of the exchange-correlation energy, potential, and kernel.All calculations were performed using a 24-processor parallel architecture.

Validation calculations
To validate the GGA ADFT kernel working equation, Eq. ( 24), we performed frequency analyses of the triplet H 2 O molecule at a not optimized geometry.To do so, we use the singlet VWN/DZVP/GEN-A2 optimized H 2 O structure, as depicted in Table 1, for our triplet frequency analyses.The obtained harmonic frequencies are listed in Table 1, too.Here, we compare analytic ADFT frequencies obtained with the LDA VWN functional and the GGA PBE functional.For both functionals, the native deMon2k functional implementation is compared with the corresponding LibXC implementation.For both implementations analytic (ANK) and finite difference (FDK), kernel results are listed.Each column displays the three harmonic frequencies of the triplet H 2 O at the optimized singlet geometry.As can be seen from Table 1, the first frequency is always imaginary (indicated by a minus sign).For the LDA VWN frequencies, excellent agreement between all methods is observed.Deviations are in the range of 1 cm -1 or below.The situation changes for the GGA PBE frequencies.Whereas the agreement of harmonic frequencies calculated with the analytic and finite difference kernels remains excellent, significant deviations between the native and LibXC kernel calculations are found.These differences are largest, up to 20 cm -1 , for the imaginary frequency and reduce with increasing frequencies as Table 1 shows.We Table 1 Analytic unrestricted ADFT harmonic frequencies (cm -1 ) of H 2 O calculated with the LDA VWN and GGA PBE functionals employing their native deMon2k and LibXC implementations.In the analytic frequency analyses, analytic kernel (ANK) and finite difference kernel (FDK) calculations are used.The H 2 O structure for these frequency analyses is depicted on the right.See text for further details attribute this to different screenings in the native and LibXC partial derivatives of the density-weighted PBE exchangecorrelation energy density.Similar deviations are found for other GGA functionals, too.
To further investigate the differences between the native and LibXC implementations of exchange-correlation functionals, we performed frequency analyses on small openshell systems.In these unrestricted calculations, the VWN/DZVP/GEN-A2 and the PBE/DZVP/GEN-A2 level of theories were used.To be unbiased with respect to optimized structure parameters, we used structure parameters from the National Institute of Standards and Technology's Computational Chemistry Comparison and Benchmark Database (NIST) [55] for the frequency analyses.
Table 2 compares the calculated harmonic frequencies of ten selected open-shell systems with each other.The table is structured in the same way as Table 1 adding mean absolute deviations (MADs) with respect to the native deMon2k ANK results at the bottom of the table.As Table 2 shows, a general good to excellent agreement between all frequencies is observed.This is confirmed by the MADs that are always below 1 cm -1 .In particular, the difference between harmonic frequencies from the native and LibXC functional implementations is reduced well below 10 cm -1 .Largest deviations are observed for 2 BeH, both for the LDA VWN (∼ 5 cm -1 ) and GGA PBE (∼ 7 cm -1 ) functionals.Furthermore, Table 2 shows that the differences between native and LibXC implementations are now in the same size range as for analytic (ANK) and finite difference (FDK) kernel calculations (∼ 5 cm -1 for 2 PH 2 ν 1 with the LDA VWN functional).This supports our previous assumption that these differences arise from different screenings in the functional implementations.We attribute the improvement of the consistency of the frequencies in Table 2 with respect to Table 1 to the molecular geometries.Whereas the molecular structure parameters for the systems in Table 2 are close to optimized structure data, the triplet H 2 O structure used in Table 1 is far from the corresponding minimum.This is also supported by the fact that the largest deviation between native and LibXC functional implementations is observed for the imaginary frequency of the triplet H 2 O example.

Benchmark calculations
To benchmark our new GGA ADFT kernel implementation, we revisited the (Pt 3 Cu) n clusters previously reported in [11,12].To this end, we re-optimized the cluster structures and performed frequency analyses with the optimized structures employing the native and LibXC implementations of the RPBE [52] functional.We used this functional because of its use in the original report.Due to the computational demand of  these calculations, we profiled some of them.To our surprise, we found that a significant portion of the computational time was spent on the finite difference calculation of the ECPs.Table 3 lists these timings in column 2.
To overcome this computational bottleneck analytic second ECP (and MCP), derivatives were implemented in deMon2k as described in Section 2.4.The timings for the analytic ECP second derivative calculations are reported in column 3 of Table 3 and graphically compared with their finite difference counterparts in Fig. 1.Speed-ups of up to almost 100 are observed for the studied (Pt 3 Cu) n clusters!As Fig. 1 shows, the newly analytic implementation of ECP sec-Fig. 1 Comparison of the computational timings for finite difference and analytic ECP second derivative calculation of (Pt 3 Cu) n clusters, with n=1 to 11.All calculations are performed on 24 cores ond derivatives exhibits a nearly linear scaling with respect to the number of ECP centers.Note that the semi-local integrator scaling is only due to the numerical radial integrator [45] because analytic angular integration is independent of the number of pseudopotential centers in the system.
The five lowest harmonic frequencies of the optimized (Pt 3 Cu) n are listed in Table 4 for the native and LibXC functional implementations.The MADs with respect to the native deMon2k ANK frequencies are given at the end of the table.Again, the analytic (ANK) and finite difference (FDK) kernel calculations are compared for both implementations.Notably, there is an outstanding agreement between the analytic and finite difference kernel calculations as shown by the Table 4 Analytic unrestricted ADFT harmonic frequencies (cm -1 ) of (Pt 3 Cu) n clusters, with n=1 to 11, calculated with the GGA RPBE functional using its native deMon2k and LibXC implementations.In the analytic frequency analyses, analytic kernel (ANK) and finite dif-ference kernel (FDK) calculations are used.The molecular structures are re-optimizations of the structures reported in [11,12] corresponding MADs that are all below 0.5 cm -1 .Deviations rarely exceed three wavenumbers, underscoring the robustness of the deMon2k implementation.In contrast, the results obtained with the LibXC interface exhibit more deviations.While deviations are generally less than two wavenumbers (maximum MAD is 2.4 cm -1 ), some frequencies display discrepancies of up to 10 wavenumbers.Additionally, using the LibXC interface with the above-mentioned methodology, two systems revealed imaginary frequencies in the case of the FDK method, namely Pt 6 Cu 2 and Pt 27 Cu 9 .To overcome this problem, we re-optimized these two clusters with a tighter structure optimization tolerance of 10 -5 a.u. and an analytically calculated start Hessian matrix.In the subsequent frequency analyses, only positive frequencies are observed.For consistency, both ANK and FDK frequencies for Pt 6 Cu 2 and Pt 27 Cu 9 listed in Table 4 were obtained with this methodology.These results show that the here-described ADFT kernel implementations are numerically stable for delicate electronic structures typical of transition metal clusters.

Conclusions
The unrestricted GGA ADFT exchange-correlation kernel working equation for perturbation dependent auxiliary functions is derived and discussed.We validated and benchmarked this kernel implementation in the framework of frequency analyses using the native and LibXC exchangecorrelation functional implementations in deMon2k.Overall, there is a very satisfying consistency between the native and LibXC frequencies.In both cases, deviations of less than a wavenumber are typically observed when comparing analytic and finite difference kernel calculations.Between the implementations, larger deviations of 5 cm -1 or more are occasionally observed.The largest deviation (∼ 20 cm -1 ) between native and LibXC frequencies is found for the imaginary mode in the non-minimum structure of 3 H 2 O.We attribute this to the different screenings in the exchangecorrelation energy density calculations.However, for (near) minimum structures, the agreement of native and LibXC frequencies is for small open-shell molecules usually excellent.
In the case of (Pt 3 Cu) n clusters, results obtained using the LibXC interface show more significant deviations, albeit generally within 2 wavenumbers, with occasional discrepancies of up to 10 wavenumbers.In the case of Pt 6 Cu 2 and Pt 27 Cu 9 , these deviations led to imaginary frequencies in the final optimized structures with the default convergence thresholds.By tightening these thresholds and using an analytically calculated start Hessian matrix, these problems were fixed and the clusters converged with the LibXC functional implementation to minima, too.The computational bottleneck due to the finite difference second derivative ECP (and MCP) integral calculations in the (Pt 3 Cu) n clusters was resolved by the implementation of corresponding analytic second derivatives.The resulting speed-ups are large and scale linearly with the number of ECP (MCP) centers.Consequently, frequency analyses of systems with hundreds of ECP (MCP) centers are now feasible with deMon2k in very reasonable times.

Table 3
Computational timings (s) for finite difference (FD) and analytic (AN) ECP second derivative calculations of (Pt 3 Cu) n clusters, with n=1 to 11.The speed-up factor refers to frequency analyses performed on 24 cores